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Abstract A single-crystal gradient plasticity 
model is presented that includes a power-law 
type defect energy depending on the gradient 
of an equivalent plastic strain. Numerical regu¬ 
larization for the case of vanishing gradients is 
employed in the finite element discretization of 
the theory. Three exemplary choices of the defect 
energy exponent are compared in finite element 
simulations of elastic-plastic tricrystals under 
tensile loading. The influence of the power-law 
exponent is discussed related to the distribution 
of gradients and in regard to size effects. In ad¬ 
dition, an analytical solution is presented for 
the single slip case supporting the numerical re¬ 
sults. The influence of the power-law exponent is 
contrasted to the influence of the normalization 
constant. 
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1 Introduction 

Continuum theory based material models are 
widely used in order to investigate the material 
behavior of microspecimen under different load 
settings like torsion (e.g., [1]). In the micron 
regime, non-classic material behavior like size 
effects occur, see, e.g., [2]. These phenomena can 
be modeled with gradient crystal plasticity the¬ 
ories by incorporating an internal length scale, 
e.g., [3]. As discussed by [4], the need for the 
incorporation of an internal length scale is a 
consequence of the coarsening error made in the 
microstructural energies of continuum theories. 
One method of incorporating internal length 
scales is the use of plastic strain gradient ex¬ 
tended work hardening laws [5]. Commonly, how¬ 
ever, classic plasticity theories are extended with 
higher-order quantities instead. These theories 
can be further classified into work-conjugate and 
non-work-conjugate theories [6]. In the work- 
conjugate theories, the free energy is extended 
by a gradient related energy contribution. This 
defect energy considers, for example, the gradi¬ 
ents of plastic slip or plastic strain. In the non- 
work-conjugate theories, however, backstresses 
on the individual slip systems are deduced from 
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physical considerations of collective dislocation 
behavior. Furthermore, in [5], it is shown that 
these two approaches are equivalent within a 
three-dimensional multislip context. Both types 
of theories can be used to predict the same kind 
of length scale effects in the material response. 
The extension of the free energy is usually ac¬ 
complished by adding a plastic contribution to 
the standard elastic free energy, see, e.g., [7]. As 
outlined by [8], an additive split of the free en¬ 
ergy into an elastic and a plastic contribution 
(that can, for example, account phenomenologi¬ 
cally for hardening in terms of the plastic slips) 
is a good first approximation for metals. In their 
work, the internal length scale is considered by 
taking into account the self energy of the dis¬ 
locations (depending on the plastic part of the 
deformation gradient). This energy is also uti¬ 
lized in [9], where it is shown that the occurring 
problem of non-differentiability in the origin of 
the domain of this energy can be circumvented 
by posing a dual formulation of the variational 
problem. 

For the incorporation of the defect energy, sev¬ 
eral approaches exist. The shared aim of these 
is the modeling of lattice distortion due to geo¬ 
metrically necessary dislocations (GNDs). It is, 
however, also outlined in [4] that the common 
association of the gradients to GNDs via the dis¬ 
location density tensor [10] is not unique. In [7], 
a geometrically non-linear theory is proposed 
using a defect energy being “for convenience” 
quadratic in the plastic part of the deformation 
gradient. A quadratic defect energy is also used 
in [11], formulated to depend on the geometri¬ 
cal dislocation density tensor within a geomet¬ 
rically linear theory. This formulation reduces 
to a quadratic form in the Burgers vector for 
strict plane strain. The discrete dislocation dy¬ 
namics results of [12] include a comparison of 
several energy formulations. Their results moti¬ 
vated [13] to reformulate the defect energy with 
dependence on the densities of screw and edge 
dislocations rather than on the Burgers tensor. 


It is shown there, that this defect energy can be 
recast in terms of the gradient of plastic slips. 
The discussion of scaling regimes by [14] includes 
this defect energy by [13] in their intermediate 
scaling regime. 

A defect energy that is quadratic in the geo¬ 
metric dislocation density tensor has also been 
used in the recent variational implementations 
of finite gradient plasticity by [15-17]. In [15], 
constitutive rate-type and algorithmic incremen¬ 
tal potentials are defined in order to formulate 
a general framework of inelasticity, applicable, 
e.g., to construct single-crystal gradient plas¬ 
ticity. In-line with this framework, a computa¬ 
tional approach to gradient plasticity (of von 
Mises-type in the logarithmic strain space), us¬ 
ing mixed variational principles, is proposed in 
[16]. This includes separation into variables of 
long-range and short-range character, i.e., of 
macro-motion and micro-motion, respectively. 
In the variational approach [17], such a separa¬ 
tion is performed for multiplicative plasticity and 
a viscous regularization technique is proposed 
to overcome the problems of classic active-set 
search for rate-independent plasticity. In addi¬ 
tion, mixed variational principles have been ex¬ 
ploited for small strains in the previous work 
[18]. 

Computational aspects of gradient plasticity the¬ 
ories are also addressed, e.g., in the works [18-29] 
and, more recently, in [30-34]. A general method¬ 
ology to incorporate the gradients of the plastic 
state variable within a suitable numerical scheme 
has been proposed in [19]. The plastic multiplier 
is considered as an independent variable, and the 
constitutive plastic equations are satisfied only 
in a weak sense. This has been extended using a 
penalty approach to allow for a C°-continuous 
interpolation of the plastic multiplier in [20]. An 
element-free Galerkin formulation has been ex¬ 
ploited for gradient plasticity in [22]. It ensures 
higher-order continuity for the shape functions 
to interpolate the plastic multiplier. 

The reformulation of strain gradient plasticity 
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presented in [21] offers benefits for numerical so¬ 
lutions: only the primary variables and their first 
gradients enter the variational formulation. This 
is an advantage compared to earlier concepts, 
e.g., [35], where the second gradients addition¬ 
ally enter the variational formulation. 

The implementation of a gradient damage frame¬ 
work [23] includes the solution of the algorith¬ 
mic consistency condition in weak form, and an 
active-set search to identify nodes in the actively 
damaged domain. Improved numerical costs are 
obtained by the implementation of gradient plas¬ 
ticity in [24], compared to the earlier works of, 
e.g., [19, 20]. In this approach, the additional 
boundary conditions need only be applied on the 
external boundaries. Thus, the increased effort 
for enforcing additional conditions on the (inter¬ 
nal) elastic-plastic boundary has been overcome. 
Discontinuous Galerkin formulations are ex¬ 
plored, e.g., in [25, 26], and an extension to finite 
strains can be found in [27]. In the formulation 
[29], it has been shown that the return mapping 
algorithm for classic plasticity can also be ap¬ 
plied to non-local plasticity such as single-crystal 
gradient plasticity, in an element-wise manner. 
A numerical treatment of a strain gradient plas¬ 
ticity theory has been presented in [30]. The 
model is of viscous type, and considers both dis¬ 
sipative and energetic contributions. This is in 
contrast to frameworks, such as the one treated 
in [29], which consider purely energetic contri¬ 
butions. Within the implementation [31] of the 
gradient plasticity theory [36], all plastic slips 
are considered as additional degrees of freedom 
and interface elements are used for the grain 
boundary discretization. Results are presented 
for the two-dimensional case. In [32], however, 
an implementation for three dimensions is dis¬ 
cussed and numerical examples are shown for 
both flow rules proposed by [36] as well as for 
an additionally proposed form. 

The approach [33] suggests using a two-step solu¬ 
tion scheme for gradient plasticity, separating the 
computations into a global solution and a local 


approximation that facilitates meshfree methods. 
It is discussed that this two-step procedure sim¬ 
plifies the implementation of gradient theories, 
with an emphasis on the enforcement of bound¬ 
ary conditions. 

Recently, a variational framework with thermo¬ 
mechanical coupling for finite strains has been 
proposed [37]. The approach [38] accounts for 
thermal annealing, and a computational study 
of this model, including viscous regularization, 
is presented in [39]. Thermal effects, however, 
are neglected in the present contribution as the 
focus is on the (mechanical part of) the defect 
energy. 

Apart from the commonly used quadratic defect 
energy formulations, non-quadratic forms have 
also been proposed in the literature. In [40], a 
defect energy is introduced that is linear in the 
accumulated GND-densities, thereby leading to 
a constant higher-order (micro) stress. A defect 
energy of more general type is used by [41]. In 
their work, several special cases are discussed 
where they distinguish between recoverable and 
nonrecoverable defect energies. Recoverable de¬ 
fect energies in this regard means that loading 
and reverse loading starting from a set of slip 
gradients leads to the same value of the defect 
energy. Nonrecoverable defect energies means, 
however, that the same value of the defect energy 
is not necessarily obtained after such a loading 
cycle. In [42], it is shown that the defect energy 
proposed by [43] (see also [44]), is recoverable but 
non-differentiable for vanishing slip gradients. It 
is shown there, that the defect energy by [40], 
which is linear in the accumulated dislocation 
densities is, however, recoverable. In [43], it is 
pointed out that the higher-order (micro) stress 
resulting from their defect energy has the form 
of a step function. This behavior is in contrast 
to the resulting behavior from the quadratic 
defect energies which lead to gradient stresses 
changing linearly in dependence on, e.g., the 
slip gradients. Furthermore, these two different 
defect energy approaches are interpreted phys- 
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ically in [43]. The dislocation self-energy leads 
to a constant gradient stress. The dislocation 
interaction-energy, however, leads to a gradient 
stress changing with slip gradients, and, thereby, 
accounting for the strain hardening induced by 
dislocation pile-ups. In [45], non-local effects are 
accounted for in the deformation of micropillars 
by considering the energy of dislocation surface 
steps, in addition to the self-energy of disloca¬ 
tions. The thermo-mechanical framework of [38] 
also uses the linear defect energy by [43] in the 
mechanical contribution to the free energy. 

In [46], linear as well as quadratic defect energies 
are reviewed, and it is shown that the latter ones 
lead to physically realistic slip profiles but also 
to unusual scaling laws. Instead, a logarithmic 
defect energy is proposed and connected to the 
statistical theory of [47]. Such a logarithmic de¬ 
fect energy is further investigated in [48] and 
compared to a linear defect energy by analytical 
calculations. In addition, numerical calculations 
with a quadratic regularization of both energy 
approaches are performed to resolve the issues 
of differentiability for vanishing GND-densities. 
A non-convex strain gradient plasticity model 
for patterning is proposed in [49]. Although a 
quadratic defect energy is used, an additional 
polynomial in terms of plastic slips is employed 
in the free energy. The choice of the value of the 
internal length scale is shown to affect the distri¬ 
butions of plastic slip. This model is extended 
in [50] using a non-convex latent hardening for¬ 
mulation by [8]. 

A generalized power-law type defect energy is 
proposed by [51] and [52]. It is formulated in 
dependence of Nye’s dislocation density tensor. 
This dependence has originally been proposed 
by [53]. The theory of [54] considers short-range 
and long-range dislocation interactions in the 
defect energy. With reference to [51], it is dis¬ 
cussed that the power exponent of the proposed 
defect energy governs the non-linearity of the 
defect energy approach. In the recent work of 
[55], it is proposed to consider the defect energy 


as a function of two invariants of Nye’s tensor. 
A power-law defect energy is proposed for this 
dependence and investigated for the choice of 
a quadratic defect energy exponent. It is also 
shown there that the experimental results by 
[56] can be fitted with a regularized logarithmic 
defect energy. 

In other recent works, it is discussed that further 
phenomena should be considered in the devel¬ 
opment of appropriate crystal plasticity models 
accounting for length scale effects. These include, 
for example, dislocation transport [57] and va¬ 
cancy diffusion in addition to coupled dislocation 
glide-climb mechanisms [58]. 

The work at hand aims at fully three-dimensional 
simulations with manageable computational 
costs. Thus, an equivalent plastic strain is used 
instead of considering all plastic slips as degrees 
of freedom. The theory is formulated for single¬ 
crystals using a micromorphic approach [59] to 
include the micromorphic counterpart of the 
equivalent plastic strain as additional degree 
of freedom in the formulation. For a brief dis¬ 
cussion of the micromorphic modeling approach 
see also [3]. The motivation to use a micromor¬ 
phic quantity is the beneficiary numerical im¬ 
plementation of this approach within the used 
framework. Gradient theories can be shown to 
be related to micromorphic theories since they 
are special cases of micromorphic theories [60]. 
This connection is constituted by constraining 
the micromorphic variable to be equal to its 
(macro) counterpart, i.e., in the work at hand, 
the equivalent plastic strain. 

The present contribution focuses on the gen¬ 
eralization of a defect energy depending on the 
gradient of this micromorphic variable. This gra¬ 
dient may be interpreted as an approximative 
measure for the densities of GNDs [61]. The 
approach of considering one equivalent plastic 
strain can, however, not be expected to give 
the same physical richness of models taking into 
account all plastic slips individually, see, e.g., 
[36]. Nevertheless, in comparison to experiments, 
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promising results have been obtained with the 
use of this quantity for single-crystals [ 1 ] as well 
as for ensembles of crystals, e.g., oligo-crystals 
[3]. The contribution at hand is considered being 
a step towards a further development of this nu¬ 
merically efficient gradient-enhanced continuum 
single-crystal plasticity model. 

It should be noted that in the present work, it is 
focused on the defect energy approach and, thus, 
on the continuum modeling of GNDs. Neverthe¬ 
less, the framework at hand would also allow one, 
with manageable effort, to phenomenologically 
include further dislocation related phenomena 
like, e.g., hardening effects stemming from sta¬ 
tistically stored dislocations. 

Outline. At first, the theoretical framework is 
presented inclnding a power-law defect energy ap¬ 
proach. Subsequently, the numerical implemen¬ 
tation within a finite element setting is briefly 
addressed. The need for a regularization method 
in the case of vanishing gradients of the mi- 
cromporphic variable is discussed and a numer¬ 
ical approach to resolve this issue is presented. 
Three-dimensional multislip finite element simu¬ 
lations of an elastic-plastic tricrystal are shown 
to illustrate the effect of the choice of the power- 
law exponent on the model behavior. Finally, an 
analytical single slip solution is presented and 
compared to the numerically obtained results. 

Notation. A direct tensor notation is preferred 
throughout the text. Vectors and 2nd-order ten¬ 
sors are denoted by bold letters, e.g., by a or A. 
A linear mapping of 2nd-order tensors by a 4th- 
order tensor is written as A = C[JB]. The scalar 
product and the dyadic product are denoted, e.g., 
hy A - B and B, respectively. The compo¬ 
sition of two 2 nd-order tensors is formulated 
with AB. The 2nd-order unity tensor is denoted 
by I. Matrices are denoted by a hat, e.g., by e. 


2 Strain gradient plasticity model with a 
power-law defect energy approach 

2.1 Motivation 

The quality of the gradient plasticity modeling 
of size effect related phenomena is significantly 
influenced by the choice of the defect energy. 
Depending on the defect energy exponent [51], it 
has been observed that the magnitude of size ef¬ 
fects of the material is varying. Common choices 
for the defect energy exponent are one and two. 
In [51], however, the choice of a linear defect 
energy is argued against. This is based on the 
argument that a linear defect energy leads to 
constant, and thus, pile-up independent gradi¬ 
ent stresses. On the path towards rehning the 
single-crystal gradient plasticity model of [3], a 
power-law defect energy is realized in the present 
contribution. This (in general) non-linear energy 
is formulated in dependence of the gradient VC 
of the micromorphic counterpart C to an equiva¬ 
lent plastic strain 7 eq. The approach taken allows 
for investigations with different exponents of the 
defect energy. 


2.2 Basic Assumptions 

In a geometrically linear framework, the po¬ 
sition of the material points of a body B 
is denoted by x. A Cartesian coordinate 
system is considered with the basis vectors 
reading { 61 , 62 , 63 }. The displacement gra¬ 
dient is given by grad{u) = dui/dxjCi ^ ej, 
and the infinitesimal strain tensor reads 
e = sym(grad (it)). The plastic strain tensor re¬ 
sults from a superposition of slips 

£P = ^A„M*, (1) 

a. 


where the symmetric part of the Schmid ten¬ 
sor is given by = sym(dQ, ® Ua). Slip direc¬ 
tions are denoted by d^, and n^, are the slip 



6 


E. Bayerschen, T. Bohlke 


plane normals of the slip systems a. For the face- 
centered cubic (FCC) crystal lattice case consid¬ 
ered here, the two possible slip directions of each 
slip system are considered separately. Therefore, 
N = 24, and the positive, non-decreasing quan¬ 
tities Aq, a = 1,..., 24 are addressed by plastic 
slip parameters, in the following. It is remarked 
however, that a treatment with 12 direction- 
ally dependent plastic slips would also be suit¬ 
able. Performing the common additive decompo¬ 
sition of the strain tensor gives the elastic strain 
tensor e® = e — . An equivalent plastic strain 

measure is introduced as 

7eq(A) = ^ / Aadf = ^ Aa, (2) 

a OL 

see, e.g., [61]. 

2.3 Principle of virtual power and field 
equations 

The principle of virtual power is used for the 
derivation of the field equations. A micromor- 
phic variable C is used as an additional degree 
of freedom in the framework (see also the intro¬ 
duction of the present work). Body forces are 
neglected, in the following. The virtual internal 
power is assumed to be given by 

5Vmt = J ^cr • de -I- TrdC + 4 du, (3) 

B 

with the virtual strain rates 6s{6u), and the 
virtual rate of the micromorphic variable It 
is remarked that and C are a priori chosen 
to be independent, i.e., = 0 (cf. [62]), and 

that the micromorphic variable C is penalized to 
be equal to its (macro) counterpart, the equiv¬ 
alent plastic strain 7 eq, see Eq. (5) and Eq. (8). 
The Cauchy stress tensor is denoted by cr, the 
stress associated to the rate C is denoted by tt, 
and the gradient stress | is associated to the 
rate Vi). All three occurring stress quantities, cr, 
TT, and are work conjugate to their associated 


kinematic quantities e, C and V^, respectively. 
It is assumed that C is a continuous quantity. 
Consequently, possible jumps of C, e.g., at grain 
boundaries in the body B, are neglected. This 
assumption is not expected to be transferable 
to plastic slips of individual slip systems, in gen¬ 
eral. In the context of an overall description of 
plastic slip, it is remarked that the calculations 
of an effective plastic strain in the experimental 
work of [63] leaded to continuous distributions 
of this overall quantity across many GBs. Thus, 
it appears to be reasonable to consider C as a 
continuous quantity in a first approach. 
Furthermore, the virtual power of external forces 
is assumed to have the contributions 

^’Pext = J i-Su da + J S'(5da, (4) 

dBt dBs 

where power is expended at the external bound¬ 
aries dBt U dBs = dB by the tractions i, and 
the microtractions S, respectively. In Eq. (4), 
it has been exploited, already, that the virtual 
rates vanish for given {m, C} at the 

Dirichlet boundaries dB^- By using the prin¬ 
ciple of virtual power, i.e., by letting the virtual 
power of the internal forces be equal the vir¬ 
tual power of the external forces, SVint = dVext, 
and applying Gauss’ theorem, the field equations 
and the Neumann boundary conditions can be 
derived (cf. appendix A). They are shown in 
Box 1. The classic quasistatic balance of linear 
momentum is supplemented by an additional 
microforce balance taking into account the gra¬ 
dient stress arising from the gradient extension 
of the free energy. Consequently, in addition to 
the Neumann boundary conditions (BCs) for 
the Cauchy stress, Neumann BCs for the gradi¬ 
ent stress are necessary, as well. It is remarked 
that contrary to the theories using all plastic 
slips as additional degrees of freedom (DOF), 
in the theory at hand, only one scalar, the mi¬ 
cromorphic variable C (penalized to be equal to 
the equivalent plastic strain 7eq), is employed 
as additional DOF. Consequently, the gradient 
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Box 1: Field equations and boundary conditions. 


Linear momentum balance 

0 

= div (cr) 

Vx e B 

Microforce balance 

TT 

= div($) 

Vx e B 

Neumann BCs for: Cauchy stress 

an 

= t 

on dBt 

Gradient stress 

$ • n 

— S 

on dB s 


stress associated to the gradient of the mi- 
cromorphic counterpart of the equivalent plastic 
strain, is identical on all slip systems. In theories 
considering all plastic slips as individual degrees 
of freedom, however, different gradient stresses 
are obtained on the individual slip systems. It 
should be noted that microhard conditions are 
employed on the interfaces between the crystals 
(i.e., on the grain boundaries) and on the exter¬ 
nal boundaries in the lateral direction (i.e., the 
two planar surfaces with minimum and, respec¬ 
tively, maximum x-coordinate). This means that 
the equivalent plastic strain remains zero, there. 

2.4 Constitutive assumptions 
2.4-1 Free energy density 

The free energy density is assumed to have the 
form 

W{e, A, C, VC) = IFe(£, £P(A)) + ITg(VC) 

+IT^(C-7eq(A)), 

with a quadratic elastic contribution 

IFe(e,eP(A)) = i(e - eP) • C[e - eP], (6) 

where C denotes the elastic stiffness tensor, and 
a generally non-quadratic gradient-related defect 
energy contribution 

VFg(VC) = IFo . (7) 

Here, Wq is an initial defect energy. For the 
case of a defect energy being linear in the dis¬ 
location density, this value is commonly taken 


to be IFo = dG, cf. [43, 64]. The shear modu¬ 
lus is denoted by G, and d ss 0.5 ... 1 is often 
assumed to be a constant. The quantity qq de¬ 
notes a normalization constant and can, thus, 
be related to the internal length scale that is 
introduced by the defect energy in the theory. In 
the case of a linear defect energy, this constant 
could be taken to be the inverse of the Burgers 
vector length, i.e., 50 = 1/^ (see, e.g., [43]). It 
should be noted that the gradient extension is 
performed here in terms of the micromorphic 
held variable C instead of 7 eq, cf. Eq. (5). The 
coupling of both quantities is accomplished with 
the energy density 

W^x(C-7eq(A)) = ^(C-7eq)^ (8) 

which penalizes deviations of C, from 7 eq. By use 
of a large penalty factor Fl^ it is achieved that 
c « 7eq (see also [3]). 

2 . 4.2 Dissipation 

The total dissipation is the difference between 
the power of the external forces and the rate of 
the free energy in the bulk, 

Aot = Jvdv = iPext - Jwdv>0, (9) 
B B 

if thermal and grain boundary ef¬ 
fects are neglected. By letting 
■Pext = Vint = /g (^(T • e + TtC + ^ • VC) dv, 
and substituting Eq. (5) in Eq. (9), the 
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dissipation T) can be reformulated as 


V = 

+ 



\ . 

3hEe .p 

dW^ . 

de y 

).s- 

de'P 

57eq 

dW^\ 


" 5114 A 

• VC > 0. 

dC ) 

5Vc) 


( 10 ) 


Introducing the abbreviation p = dW^/djeq = 
—dW^/d(, assuming the stresses cr, tt 
and ^ to be purely energetic, and substitut¬ 
ing dWe/de^ = —cr, leads to the reduced 
dissipation inequality 


Here, 70 is the reference shear rate, is the drag 
stress, and p is the rate sensitivity exponent. The 
initial yield stress of the slip systems is denoted 
by Tq. It is remarked that on all slip systems a, 
the term div (^) is acting identically. This is due 
to the usage of the gradient of the micromorphic 
variable VC in the defect energy Eq. (7), instead 
of considering the gradients of the slips of all 
slip systems individually. 


3 Finite Element Implementation 


= (T • eP -P7eq > 0. (11) 

Substituting Eqs. (1) and (2) in Eq. (11) gives 

V = J2iTc.-p)K>0, ( 12 ) 

cx 

with the resolved shear stresses Ta = cr ■ AT®. 
Eurthermore, based on [11], it is assumed that 
the bulk dissipation is induced by the slip sys¬ 
tem contributions, i.e., by the dissipative shear 
stresses and the slip parameter rates Aq,, 

v = Y,TtK. (13) 

oc 

Consequently, the dissipative shear stresses can 
be expressed by 

Ta=Ta-p. (14) 

Einally, using p = —tt from Eq. (10) and the 
microforce balance (tt = div (C), see Box 1), gives 


Ta=Ta+ div (C) . (15) 

2.4-3 Flow rule 


The flow rule is assumed to be of an overstress 
type, formulated in the plastic slip parameter 
rates A^, 


Aa = 7o 


T'a-'To 


C \ P 


. / Ta -b div (C) - 
= 7o 


c \ P 


(16) 


3.1 Linearization of the variational form 

The theory is implemented with finite elements 
(FE) in an in-house EE-code. An enhanced time- 
integration algorithm is used that allows for 
large time-steps compared to numerical standard 
schemes (see [65]). In order to discuss the key 
concepts in the implementation of the power-law 
defect energy approach, Eq. (7), the lineariza¬ 
tion of the principle of virtual power is addressed, 
briefly. The linearization of Eqs. (3) and (4) in 6 u 
yields 


B 


dcr , ^ , da ^ , 

+ •^edu = 


— J a-5edv-\- J t-Suda 'iSu. 

B dBt 


(17) 


Eurthermore, the linearization of Eq. (3), and 
Eq. (4) in SC,, considering additionally tt = —p, 
gives 


B 

+ ^[V(AC)]-V(dC)) du = 


J (^-pSC + C-^{SO^dv+ J SSCda- 

B dBn 


(18) 



Power-law Defect Energy in a Single-Crystal Gradient Plasticity Framework 


9 


The numerical integration of Eq. (17) and 
Eq. (18) is outlined in [65] in a more general 
formulation, including grain boundary contribu¬ 
tions. 

Obviously, the EE implementation of Eq. (17) 
and Eq. (18) requires the computation of the gra¬ 
dient stress 4 = 91Tg(VC)/9VC and its deriva¬ 
tive Using the power-law defect energy 

formulation, Eq. (7), these quantities can be 
expressed by 


dW^jS/Q 

avc 


”*]Sivcr-" 

yo 


VC, 


(19) 


3.2 Regularization of the power-law defect 
energy approach 

A standard regularization approach is under¬ 
taken here, e.g., [66]. Regularized defect energy 
associated quantities are denoted by a (•), in 
the following. The regularized defect energy is 
introduced as 


lUg(VC) = lUo 




.2\2 


90 


9o 


(23) 


+ m(m-2)^lVCl'”-"VCOVC, 

9o 


( 20 ) 


respectively. From Eq. (19) and Eq. (20) it be¬ 
comes apparent that for an arbitrary power ex¬ 
ponent TO, the gradient stress in this formulation 
is neither defined nor differentiable in the case of 
vanishing gradients VC, i.e., for jVC] —> 0. How¬ 
ever, for the special case of a quadratic defect 
energy with to = 2 , the gradient stress reads 

4 = 2^VC, (21) 

9o 

and its derivative with respect to VC is 


aVC ' 


( 22 ) 


Consequently, for this special case, the gradi¬ 
ent stress is defined and differentiable, even 
for vanishing gradients VC. It is remarked that 
the quadratic defect energy leads to a gradient 
stress ^ that is linear in the gradient VC. Choos¬ 
ing TO = 1, however, yields a constant gradient 
stress 

In order to use the power-law defect energy with 
varying exponents to, a regularization of the 
power-law defect energy, Eq. (7), for vanishing 
gradients VC, is outlined in the following. 


The numerical parameter e is chosen sufficiently 
small but non-zero (see Section 4.2), and leads 
to finite values of the derivatives of Eq. (23), 
even in the case of vanishing gradients VC. 

The gradient stress corresponding to the em¬ 
ployed regularization is 


€ = 


9IUg(VC) „_IUo 


avc 


= TO—(IVOR¬ 
S'™ 


, ^-1 


VC, 


(24) 


and its derivative with respect to the gradient is 
given by 


+ to(to - 2)^(1VC1" + e")^”"vc 0 VC. 

9o 


(25) 


For the special case to = 2, this regularized 
formulation recovers the non-regularized 
defect energy, i.e., C = 2HoVC/5o, 
d^/dVC = 2WoI/gQ, respectively. If the reg¬ 
ularization parameter e is set to vanish, the 
original formulation, Eq. (7), of the power-law 
defect energy is recovered, i.e., C = and 
9C/5VC = 9C/5VC, respectively. 
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3.3 Numerical time integration and algorithmic 
tangents 

The numerical time integration requires the de¬ 
termination of the stresses in Eqs. (17) and (18) 
and the computation of the algorithmic tangents. 
These computations follow along the scheme 
discussed in [3], but the grain boundary con¬ 
tributions are neglected since only microhard 
conditions are considered on the grain bound¬ 
aries throughout the present work. Furthermore, 
isotropic hardening of Voce-type is not consid¬ 
ered in the work at hand, and thus, p does not 
have an isotropic hardening contribution /? as 
in the above mentioned work. An implicit Euler 
scheme is applied to Eq. (16). For details on this 
it is referred to [65]. 

4 Numerical results 

4.1 Material model, geometry, and boundary 
conditions 

The material considered throughout the finite 
element simulations in this work is, for simplic¬ 
ity, chosen to be elastically isotropic with ma¬ 
terial constants similar to aluminum, see Ta¬ 
ble 1. Fully anisotropic plasticity, however, is 


Table 1: Model parameters for tensile test simu¬ 
lations of aluminum-like tricrystals. 


G 

V 

Wo 

CK 

p 

25 GPa 

0.3 

aO 

0.5 

20 

70 



e 

Hyc 

10“^ 1/s 

33.5 MPa 

1 MPa 

10“® 

10® MPa 


accounted for in a multislip setting for FCC ma¬ 
terials. A tricrystal composed of cubic grains 
(see Fig. 1) is simulated to demonstrate the be¬ 
havior of the model for different defect energy 
exponents m. The same initial defect energy Wq 
is used in all simulations. At first, the three 


grains are considered to be in a standard crystal 
lattice orientation (i.e., they are oriented (100) 
with respect to the cc-axis of the Cartesian x- 
y-z-system). In order to allow an assessment 



Fig. 1: Spatial discretization of tricrystal with 
elastic bounding grains (blue) and line segment 
along the central axis (white). 

of the power-law defect energy approach, the 
behavior of the occurring gradients of the micro- 
morphic variable C close to the grain boundaries 
is isolated by setting the two bounding grains 
to behave purely elastic. Subsequently, a fully 
elastic-plastic tricrystal is discussed, addition¬ 
ally, where the two grain boundaries do not allow 
any dislocation transmission across, i.e., they are 
treated as microhard boundaries. In the theory 
at hand, these microhard boundaries can be in¬ 
cluded as a subset of an arbitrarily chosen grain. 
No additional grain boundary terms need, there¬ 
fore, to be considered in the theoretical model 
(cf. also the discussion on this choice in [65]). 
The tricrystal of dimensions 3x3x9 pm 
(see Fig. 1) is loaded strain controlled up to 
an overall strain of e = 0.005 under tensile load¬ 
ing along the a;-axis. Lateral contraction is al¬ 
lowed for by the boundary conditions on the two 
bounding planes at a; = 0 pm, and x = 9 pm. In 
the following, distributions of the micromorphic 
field variable C are investigated for the differ¬ 
ent elast-plastic cases under consideration. The 
line segment highlighted in Fig. 1 is, therefore, 
utilized to obtain plots of the field distribution. 
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4.2 Convergence of results 

The convergence of the numerical results is dis¬ 
cussed twofold, here. Besides convergence in re- 



Fig. 2: Final nominal stress at f = 1.0 s vs de¬ 
grees of freedom for tricrystal with elastic bound¬ 
ing grains, different spatial discretizations, and 
different numerical regularization parameters e. 


gard to the spatial discretization of the finite 
element mesh, the convergence with respect to 
the regularization parameter e is investigated. 
The defect energy exponent used in this section 
is TO = 1.1, and the normalization constant is 
go = 450.21/|i.m. As can be seen in Fig. 2, the 
regularization parameter e is of negligible influ¬ 
ence on the final nominal stress at the end of the 
loading. It also does not influence the conver¬ 
gence behavior regarding the spatial discretiza¬ 
tion. In Fig. 3, it is shown that a sufficiently small 
choice of the regularization parameter leads to 
coinciding distributions of C. along the x-axis of 
the central grain for a further decrease in e. 



Fig. 3: Distribution of micromorphic variable C in 
central grain of tricrystal with elastic bounding 
grains for discretization with 349804 DOF, and 
different numerical regularization parameters e. 

4.3 Results for the defect energy with different 
exponents m 

4-3.1 Tricrystal with elastic bounding grains 

In order to investigate the model behavior for 
different exponents to in the defect energy ap¬ 
proach, Eq. (23), two approaches are undertaken. 
At first, numerical results are obtained with a 
subsequent comparison to results from an an¬ 
alytical solution, see Section 5. For a first in¬ 
vestigation, three choices of the exponent to 
in the defect energy approach are considered. 
For each choice, a corresponding normalization 
constant go is determined by simulations such 
that the final nominal stress value at the final 
deformation strain is identical for the choices 
of TO (and go, respectively). Consequently, for 
all considered exponents, the plasticity carried 
by the simulation specimen is identical after the 
final time step. The choices for to and the corre¬ 
sponding normalization constants considered are 
shown in Table 2. The behavior of the model can 
then be compared for the cases of a quadratic 
defect energy, a defect energy that is closer to 
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Nominal strain e 



t = 0.48 s 
t = 0.40 s -- 
t = 0.32 s ^ 


i = 0.48 s ^ 
t = 0.40 s ^ 
t = 0.32 s 


t = 0.48 s ^ 
t = 0.40 s 
t = 0.32 s ^ 


Fig. 4: Left; Stress-strain curves for tricrystal with elastic bounding grains and different defect 
energy exponents m. Right: ^-distributions in the central grain along the a;-axis for three time steps. 


being linear, and a defect energy exponent that 
is in the range between these two cases. It is re¬ 
marked that the employment of a defect energy 
with exponent m = 1.0 would require further 
regularization due to the occurring jump in the 
gradient stresses at the beginning of the defor¬ 
mation. Instead, a value of m is investigated 
that is sufficiently higher to not necessitate ad¬ 
ditional regularization. The stress-strain curves 
that are obtained after this procedure are de¬ 
picted in Fig. 4 (left). While the quadratic de¬ 
fect energy leads to a linear hardening behavior, 
choices of lower exponents m, however, reduce 
the slope of the (in general non-linear) hard¬ 
ening observed. The field distributions of the 
micromorphic variable C along the central line 
segment in the elastic-plastic grain are depicted 
for the three choices of m in Fig. 4 (right). It 


can be seen there, that the common choice for 
the exponent of m = 2.0 yields a hyperbolic pro¬ 
file for C in the central grain. A choice closer 
to m = 1, however, leads to a more plateau-like 
distributed C, in the central grain. The gradients 
close to the grain boundary rise with decreas¬ 
ing m. The full evolution of the C-profiles is, 
furthermore, shown in Fig. 5. For the quadratic 
defect energy, a hyperbolic distribution of C de¬ 
velops. In contrast, the lower choices of m lead 
to the initiation of a region with constant () that 
proceeds throughout the evolution. The lower 
the choice of m, the more pronounced is this 
region. Due to the steeper gradients developing 
for these m, the maximum value that C, achieves 
is significantly less than for m = 2.0. The same 
“amount” of plasticity is carried by the central 


Table 2: Used defect energy exponents m and 
corresponding normalization constants go with 
respective internal length scales 1/go- 


771 

2.0 

1.5 

1.1 

9q in 


8.25 

43.5 

450.21 

l/flo 

in [Am 

1.2 X 10“^ 

2.3 X 10“^ 

2.2 X 10“^ 














Power-law Defect Energy in a Single-Crystal Gradient Plasticity Framework 


13 


CQOOC^fOlOt'COOCM 

--hOOOOCOOOO 





3 3.5 4 4.5 5 5.5 6 3 


X in |i.m 


4 4.5 5 

X in [im 


6 3 


4 4.5 5 

X in 


5.5 6 



Fig. 5: Evolution of micromorphic variable in the central grain of tricrystal with elastic bounding 
grains. 


grain, however, at the end of the last time step 
in each case. 

4-3.2 Size effects for tricrystal with elastic 
bounding grains 

In addition, the size effects occurring for the dif¬ 
ferent defect energy exponents are investigated 
with the aforementioned tricrystal. The same 
model parameters are considered as in the previ¬ 
ous example, the length I of the cubic grains is, 
however, varied. In Fig. 6, the different size effect 
behavior is shown for three grain sizes 1. Smaller 
grains respond stiffer for all three values of to, 
and the magnitude of the size effects is clearly 
influenced by the choice of the defect energy 
exponent and the corresponding internal length 
scales, which are also indicated in Fig. 6. It is 
remarked that the model behavior is influenced 
both by the choice of the defect energy exponent 
and the internal length scale. An increase in the 
overall yield strength becomes more pronounced 
for smaller to- values. Both phenomena are dis¬ 
cussed, in more detail, in Section 5. For constant 
defect energy exponent and length scale, the 
tricrystal responds stiffer with decreasing grain 



Nominal strain e 


Fig. 6: Stress-strain curves for tricrystal with 
elastic bounding grains for different grain sizes I 
and defect energy exponents to. 

size 1. A higher magnitude of the size effect is 
observed for the quadratic choice, compared to 
the smaller choices of the defect energy expo¬ 
nent. The smaller choices of to, however, show 
a less pronounced size effect. In addition, the 
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corresponding distributions of the micromorphic 
variable ( are plotted for the final time step of 
each simulation in Fig. 7. It can be observed, 
there, that the size of the grains affects the dis¬ 
tribution of significantly. Larger grains tend 
to show shorter plateaus of ( and larger values 
of ( in the center of the grain. The classic size 



X/Xn 


Fig. 7: ^-distributions in the central grain along 
the cc-axis plotted over normalized coordinate 
2;/a;niax at the final time step for different grain 
sizes I and defect energy exponents m. 




effect of “smaller responding stronger” can be 
observed in both Fig. 6 and Fig. 7. In detail 
this means that the smaller the grain size, the 
smaller is the ” amount “ of plasticity carried by 
the central grain, see Fig. 7. 

4-3.3 Tricrystal with elastic-plastic grains 

As an additional numerical example, a tricrys¬ 
tal consisting of cubic grains behaving all fully 
elastic-plastic is investigated. This shall demon¬ 
strate the capability of the employed defect en¬ 
ergy approach to be used in fully elastic-plastic 
three-dimensional multicrystalline simulations. 
The material parameters are chosen to be the 


Fig. 8: Field distributions of C for tricrystal with 
elastic-plastic grains (top), and C-distribution in 
the central grain along the x-axis (bottom) for 
defect energy exponents m. 

same as in the previous numerical examples. 
However, the crystal orientation of the central 
grain is altered by a rotation of 6° around the z- 
axis (see Fig. 1). The resulting held distributions 
of C are shown in Fig. 8 (top) for all three choices 
of TO. Clearly, the width of the gradient affected 
zone close to the grain boundaries at cc = 3 p.m, 
and at a: = 6 p.m, respectively, is influenced by 
the defect energy exponent to. Lower choices 
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lead to steeper gradients that are distributed 
over smaller regions. This can also be seen in 
Fig. 8 (bottom), where it can be observed, as 
well, that the different orientations of the grains 
lead to differing gradients of C,. For m = 1.1 this 
becomes especially apparent with a very sharp 
transition in the C-distribution, from the plateau 
in the middle of each bounding grain towards 
the grain boundaries. 


5 Exact Solution of a Laminate for 
Single Slip 

In order to confirm the distributions of the 
micromorphic variable C obtained with the 
power-law defect energy, as shown in Section 4, 
an analytical solution is briefly presented. By 
Eqs. (7) and (10), the gradient stress is given by 


I 


dW^jS/Q VFo VC 

avC go \ 90 ) |VC| ■ 


(26) 


In the following, single slip for simple shear is 
considered in the ideal case of coinciding micro¬ 
morphic variable and equivalent plastic strain, 
i.e., C = 7 eq = 7 , and 7 = 7 (x). Consequently, 
V 7 ~ 63 ; and the x-component of ^ reads 




ix=m 


Wq 

50 



(27) 


The term div (^) in Eq. (15) reduces to 

ax 


(28) 


which, in combination with Eq. (27) gives 


—— = (m — 1 )to^ 5 - - 

da; 5o \ 5o 


d7|\ ™ ^ q|±y 


d|SI 

da; ’ 


(29) 


where 


d| 


d7 I 
dx I 


dx 


= sg 


dyN d^7 


dx/ dxdx 


(30) 


For the stationary single slip case = Tq, 

Ta = T, and, thus, Eq. (15) yields 

^ = To® - T = -At, (31) 

where At denotes the difference between the 
resolved shear stress and the critical shear stress. 
Combining Eqs. (29) and (31) results in 


—At = {m 





(32) 


with Xn,^in ^ 0 ^ a;max and Xmin — ^max- A so 
lution for 7 (x) is obtained, for example, by con¬ 
sidering sg (dy/dx) = -|-1 for —Xmax < x < 0 in 
conjunction with the BCs d 7 /dx(x = 0) = 0, 
and 7 (x = —Xmax) = 0. The solution is 


7 (x) 


lTo(w - 1) 

At 




(33) 


Analogously, the solution for 0 < x < Xmax can 
be obtained by considering sg (dy/dx) = — 1 , 
and the BCs dy/dx(x = 0) = 0, 
and y(x = Xmax) = 0 , respectively. 

In addition to the y(x)-distributions, the 
shear stress evolution in dependence of the 
average plastic slip, y, is derived. Therefore, the 
following average is considered 



0 

J y(x)dx. 




(34) 


on the interval x G [xmin 7 0 ] with Xmin = —a;max 
and Ax = —Xmin- The evaluation of Eq. (34) 
considering an average shear stress difference 
of At = At leads to 


lFo(m - 1) 
At 


ATgo\'"-'^ m 
Worn) 2m-1^“’"’' ’ 

(35) 
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Solving Eq. (35) for At yields 


Z\t(7) = t-t^ 


Worn / 2 m— 1 _\™ ^ 
(36) 


If, additionally, the case m = 1.0 is considered, 
Eq. (36) reduces to 


C ^0 , ^0 , c 


(37) 


In order to compare shear distributions for differ¬ 
ent defect energy exponents, it is assumed that 
the same average shear 7 = 0.01 is present for 
the different defect energy exponents m. Then, 
the corresponding values of the dimensionality 
constants go for different m can be obtained 
from Eq. (35). This gives the values indicated in 
Table 3. In Fig. 9 (left), the normalized distri- 

Table 3: Different defect energy exponents m and 
corresponding normalization constants go with 
respective internal length scales 1 /go, values of 50 
obtained from Eq. (35) under the assumption 
that 7 = 0.01 for all m. 


m 

2.0 

1.5 

1.1 

Qo in 

l/nm 

0.39 

0.83 

2.07 

1/90 

in }a.m 

2.56 

1.2 

0.48 


butions of 7 (x) are depicted for different defect 
energy exponents m and assuming a shear stress 
difference of At = 5 MPa. All other material 
parameters used are identical to the parameters 
in the preceding section. It can be observed that 
the smaller the defect energy exponent m and 
the internal length scale 1 /go, the higher are 
the gradients of 7 ( 0 :) close to the boundaries. A 
quadratic defect energy leads to a hyperbolic 
distribution ^{x). For choices of m closer to one, 
however, more plateau-like distributions 7 ( 0 :) are 
obtained. The resulting evolution of the average 
shear stress difference Af is plotted for the plas¬ 
tic range in Fig. 9 (right). It can be seen, there, 
that a quadratic defect energy leads to uniform 
linear hardening, while choices of m closer to one 
give highly non-linear, at the onset of plasticity 
more pronounced, but subsequent less hardening, 
compared to m = 2 . 0 . 


For this case, the influence of the grain size /, 
and the dimensionality constant go, on the av¬ 
erage shear stress become clear. The average 
shear stress r scales inversely linear with both, 
I and go- For exponents m > 1, however, these 
dependencies are non-linear, see Eq. (36). This 
explains the increase in yield stress occurring for 
smaller values of m in conjunction with smaller 
grain sizes I that can be observed in the numeri¬ 
cal results, e.g., in Fig. 6 . 

In the following, the influence of the defect en¬ 
ergy exponent m is separated from the influence 
of the internal length scale 1 /go for the illustra¬ 
tive example above. A variation of the values 
of the dimensionality constant go from Table 3 
by ±5% is performed. In Fig. 10, the resulting 
normalized average shear distributions (left) and 
the resulting normalized average shear stress dif¬ 
ference courses (right) are depicted. The courses 
from Fig. 9 are indicated as a reference by solid 
lines. From Fig. 10 (left), it becomes evident that 
the resulting shear is highly non-linear in the 
variations of m and go, respectively. Changes in 
the defect energy exponent m further amplify 
variations in the gradients of 7 when changing 
the internal length scale 1/go- For constant m, 
an increase in 1 /go leads to more pronounced 
hardening behavior, see Fig. 10 (right), consis¬ 
tent with less plastic shear (see Fig. 10 (left)). 
Changing the defect exponent m, however, gives 
rise to substantially different shapes of the re¬ 
spective hardening courses, see Fig. 10 (right). 


6 Discussion of results 

The power-law defect energy approach presented 
has been regularized to avoid the issue of vanish¬ 
ing gradients in the mathematical formulation. It 
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Fig. 9: Normalized 7 -distribution of analytical solution (left) for single slip simple shear, plotted 
along the normalized coordinate. Identical 7 for all three values of m. Normalized average shear 
stress difference (right), plotted in the plastic range over normalized average plastic slip 7 . 



m 11 go (M-m) 

2.0 2.44 . 

2.0 2.56 - 

2.0 2.70 . 

1.5 1.15 

1.5 1.20 - 

1.5 1.27 . 

1.1 0.46 . 

1.1 0.48 - 

1.1 0.51 . 



ihrr. 


Fig. 10: Normalized 7 -distribution of analytical solution (left) for single slip simple shear, plotted 
along the normalized coordinate. The internal length scale is varied as indicated. Normalized average 
shear stress difference (right), plotted in the plastic range over normalized average plastic slip 7 . 


is shown in the present work that the regulariza¬ 
tion parameter used has to be chosen sufficiently 
small in order to be of negligible influence on the 
mechanical response in the gradient plasticity 
simulations. The choice of the defect energy ex¬ 
ponent m has several influences on the observed 
mechanical response. First of all, the resulting 


hardening stemming from the defect energy is 
influenced by the exponent choice. Furthermore, 
the gradients of the employed micromorphic vari¬ 
able close to grain boundaries are affected by 
the choice of m. Smaller choices of m lead to 
larger gradients. Such changes in the distribu¬ 
tion of plastic strain can also be achieved by 
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modifying the energetic length scale employed, 
see [54] and also [49]. The magnitude of size ef¬ 
fects is significantly affected by the defect energy 
exponent choice. Smaller choices of m lead to 
less pronounced size effects. It should be noted, 
however, that besides the energetic influences on 
the magnitude of size effects, also dissipative ef¬ 
fects have been investigated in the past, e.g., by 
[51]. For a fully elastic-plastic tricrystal ensem¬ 
ble of single-crystals, the benefit of the ability to 
use different energy exponents m in the model 
becomes especially apparent. With smaller m, 
larger gradients and, thus, larger pile-ups can be 
modeled close to (grain) boundaries, compared 
to larger values of m. The analytical solution 
presented confirms the influence of the choice 
of the defect energy exponent m on the stress- 
strain curves as well as on the distribution of 
plastic slip and its gradients, respectively (com¬ 
pare Fig. 4 and Fig. 9). In addition, the influence 
of the internal length scale has been shown exem¬ 
plary for the analytical solution. The resulting 
model behavior is influenced by both the defect 
energy exponent and the internal length scale, 
see Fig. 10. 


the work-hardening courses. On the other hand, 
changes of the normalization constant lead to 
different hardening and shear profiles for con¬ 
stant exponent but the qualitative shape of both 
the shear profiles and the hardening courses is 
maintained. The used defect energy approach 
can be utilized on the continuum scale to investi¬ 
gate the choice of the exponent in the modeling 
of phenomena like pile-ups of dislocations close 
to grain boundaries. Furthermore, it is demon¬ 
strated in the work at hand that different magni¬ 
tudes of size effects can be investigated with the 
underlying approach, depending on the defect 
energy exponent and the internal length scale. 
A regularization technique to resolve the issue 
of vanishing gradients is presented. An investiga¬ 
tion of the choice of the regularization parameter 
shows no significant dependence of the model 
response on this parameter, if it is chosen suffi¬ 
ciently small. The employed regularized defect 
energy approach seems to be promising for fu¬ 
ture investigations with the gradient plasticity 
model at hand. 

A Derivation of field equations 


7 Conclusion 

The generalized defect energy approach of the 
gradient plasticity model presented in this work 
allows one to model different types of distribu¬ 
tions of the micromorphic field variable C and 
its gradient, respectively. This can be performed 
by changing the exponent of the defect energy 
and the associated normalization constant. An 
analytical solution for single slip supports the 
numerically obtained results. It is used to show 
the influence of the defect energy exponent con¬ 
trasted to the influence of the normalization 
constant. The resulting mechanical response is 
non-linear in the normalization constant due to 
the power-law defect energy approach. On one 
hand, changes in the exponent alter the shape 
of the spatial shear profiles and the shape of 


Application of the principle of virtual 
power SVint = SVext, i.e., letting the right-hand 
side of Eq. (3) being equal to the right-hand side 
of Eq. (4), yields 


c) dv= f 

dSt 


(7 - 6e H- ^ • V^C) dt; = / t - 5uda 


j E da. 


(38) 


Then, by substituting e = sym(grad (tt)) in Eq. (38) 
and applying two forms of the divergence theorem. 


J A ■ grad {b) dv = — J div (A) • 6 dn -|- J An ■ b da, 


(39) 


and 


J a • grad (c) da = — y div (a) c da-|- J (a-n)cda, 
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(40) 

neglecting any jumps, the following is obtained after 
regrouping of terms 

div (cr) ■ Siidv + j (crn — i)-5u da 
B eBt 

-|- J (tt — div (^)) dn-|- J ■ n — E)5^ da = 0. 

B aBs 

(41) 

Requiring the left-hand side of Eq. (41) to vanish for 
arbitrary virtual rates Sii, SC, yields the field equations 
in Box 1. 
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